Decomposition of spectral density in individual eigenvalue 
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The eigenvalue densities of two random matrix ensembles, the Wigner Gaussian 
matrices and the Wishart covariant matrices, are decomposed in the contributions 
of each individual eigenvalue distribution. It is shown that the fluctuations of all 
eigenvalues, for medium matrix sizes, are described with a good precision by nearly 
' normal distributions. 
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INTRODUCTION 



Borrowing from J. Wishart the random matrices used by statisticians to construct their 
ensemble of covariant matrices [lj], E. Wigner, in the late fifties, introduced the ensemble 
of Gaussian matrices of random matrix theory (RMT) (2|. Spectral properties of these 
two ensembles are characterized by correlations generated by the repulsion among the 
levels and their properties are directly connected to two classical polynomials, the Hermite 
ones, in the Wigner case[3], and the Laguerre ones, in the Wishart case. Both ensembles, 
and specially Wigner's, have proved to be important matrix models. By the same time, 
M. Girardeau proved the mapping theorem that states that in ID the properties of a gas 
of impenetrable bosons are the same of a gas of fictitious spinless fermions |4j. In the 
limit of negligible size, the only interaction among the atoms that remains is that it is 
forbidden for two atoms to occupy the same site. This generates a repulsion among them 
in an entirely similar way of what occurs among the eigenvalues. Placed in a harmonic 
trap, it follows from this equivalence that the ground state wave function of the atomic 
system is just the joint probability distribution of RMT eigenvalues jE]? ]. Therefore, 
under these conditions, the atoms of the unidimensional boson gas constitute a physical 
realization of eigenvalues of the matrix ensemble. 

In the last decades, there has been in RMT studies a great interest in the behavior of 
individual eigenvalues at the edge of spectra. This interest followed the many applications 
that the probability distribution derived by C. Tracy and H. Widom which describes the 
behavior of eigenvalues at the border of the Gaussian ensemble 0] have found. The 
distributions of the largest and the smallest eigenvalues of Wishart matrices have also 
been derived [7|. A salient feature of these edge distributions is their asymmetry that can 
be understood as an effect of the unbalanced repulsion that eigenvalues at the border 
are submitted from the other eigenvalues. Here we are interested in investigating the 
behavior of eigenvalues not only at the edge but also at the bulk of the spectrum. This 
question has been aroused some years ago in connection with the behavior of eigenvalues 
of the two-body random ensemble as compared with the RMT ones[8j]. In the boson gas 
terms, this decomposition is equivalent to determine how individual atoms behave. 
First we remark that, in general, the density of a set of random variables can be 



expressed as a sum of the individual distributions of each variable considered in an ordered 
sequence. Of course, this just translates the fact that a variable found at a given position 
has to be or the first, or the second, or the third, and so on, of the ordered sequence. 
Therefore, if the distribution of each variable is determined, a decomposition of the 
density follows. Intuitively, it is reasonable to expect that, at the bulk of the spectrum, 
if not exactly at least approximately in the case of large matrices, eigenvalue fluctuations 
will be normally distributed and, in fact this has been proved in the limit when the 
size of the matrices goes to infinity 0, 10]. As a consequence, the exact expression of the 
eigenvalue density should allow a decomposition in terms of nearly Gaussian distributions. 
It is our purpose to show that this decomposition exists and to determine its parameters. 
This will complement the exact results obtained in the asymptotic limit when matrix 
sizes go to infinity To do this, we note that the exact density of the eigenvalues of 
these two ensembles show tiny fluctuations around the average density. These wiggles 
correspond to the peaks of the individual eigenvalue distributions, that is to the average 
positions of the eigenvalues. To find these positions we separate in the eigenvalue density, 
a smooth leading order term of the fluctuating term which vanishes in the limit of large 
matrix size. The locations of the wiggles are then determined in the fluctuating term. 
Comparing then this term with the one of the same order in the asymptotic of the 
Gaussian decomposition, the dependence of the Gaussian variances on their positions 
along the spectrum is determined. 

The Wigner and the Wishart ensembles are characterized by Dyson index (3 that 
takes the values 1, 2, 4 for the three symmetry classes,respectively, the orthogonal (GOE) 
of real symmetric matrices, the unitary (GUE) of complex hermitian matrices and the 
symplectic (GSE) whose elements are real quaternion numbers. A /3-ensemble has been 
proposed that generalizes, for arbitrary positive values of (3, the Wigner and the Wishart 
ensembles [HI]. In particular, in the limit ft — > oo, the spectrum gets frozen with eigen- 
values located at the position of the zeros of the Hermite (Laguerre, in Wishart case) 
polynomials. As (3 decreases, the eigenvalues start to vibrate around those positions. 
Using perturbation theory, a Gaussian decomposition of the eigenvalue density for large 
but finite (3 was derived for both ensembles |l3|]. Here we are extending this decomposi- 
tion for small values of (3. For large individual distributions do not overlap while in our 
case there is overlapping between neighboring distributions. In the opposite limit (3 — > 0, 
the eigenvalues become a set of uncorrelated variables and, in this case, there is strong 
overlapping among the individual distributions. 



II. UNCORRELATED VARIABLES 

Starting with the case of uncorrelated variables, we consider a set of N independent 
and identically distributed (i.i.d.) variables Xi with i — 1,2, N uniformly distributed 
in the interval (—N/2, N/2). By simple probabilistic argument, the density probability 
of finding a variable at t with n others greater and the N — 1 — n others smaller than it, 
is 

1 ' ' n\(N-n-l)>. \2 N) \2 + N, 
With n = 0, 1, , ...,N — 1, this family of beta distributions gives an exact description of the 
order statistics of the set of variables in which the nth function, F(n, t), corresponds to the 



density distribution of the (n+l)th largest variable. Immediately, we verify that summing 
up all of them, the unit density is recovered as it should, namely p(t) = ^F(n, t) = 1. 
From ([!]), it follows that these distributions have first moment and variance given by 



iV(iV - 1 - 2n) 
2(N + 1) 



(2) 



and 
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(n+ l)(N-n)N 2 
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(N + 2)(N + 1) 2 ' 
respectively 

By taking the limit N — > oo keeping n fixed, the order statistics of the largest variables 
is obtained. In this case it is more convenient to express the distributions using as variable 
the distance y = t — N/2, to the right border in terms of which, the distributions F(n, y) 
converge, for large N, to 



F(n,t) 



_ y )n 



exp(y). 



(4) 



This set of functions are known to give the density distributions of the largest variables 
of an i.i.d. sequence with a compact support [141]. In particular, the first one, F(0,y) = 



exp(y), is the Weibull distribution [14j, |15|. From 



JIJ), we find that on the average the 
n + 1). 



(n + l)th variable is located at the position y = 

Consider now the situation in which, in the same limit of a large number of variables 
the ratio is kept fixed. In this case, the beta distributions converge to the Gaussian 



F(n,t) 
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exp 
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(5) 



that give the order statistics of the variables at the bulk of the sequence. Taking the 
limit of large N in the expressions ([2]) and ([3]) the first moment and the variance become 
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and 
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The above equation shows that the variance of the individual distributions at the bulk 
scales as a ~ \^N. The case of i.i.d. variables with an arbitrary symmetric density dis- 
tribution p(x), can be mapped into the above case by the transformation t = dx'p(x'). 
The density distribution of the (n + l)th variable of the ordered sequence is then 



F(n, x) = p(x) 



(N-l)\ 



n\(N 



n 



1)! 



t(x) 
N 



t(x) 
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~1 N-n-l 



In Fig. 1, it is shown the individual contribution decomposition for a set of iV 
variables sorted from the common Gaussian distribution 



(8) 
20 i.i.d. 



p{x) = —= exp(— x 2 ). (9) 
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In this case 

t(x) = yerf(x), (10) 

where erf(x) is the error function. 

It is seen that there is a strong overlapping among the individual distributions that 
reflects the dependence of the variances with the square root of the number N of vari- 
ables. The largest variable is expected to be distributed according to the Gumbel density 
distribution jlq 

F G umbei(z) = exp[- exp(-z)] exp(-z) (11) 

in the new variable 

N 

z = -ln[y erfc(x)], (12) 

where erfc(x) is the complementary error function. Indeed, one can see that this is true 
from the very good agreement exhibited in Fig. 2, in which the largest variable in a 
sequence of iV = 100 Gaussian variables is compared against the above Gumbel density 
distribution, Eq. (ITT)) . 



III. CORRELATED EIGENVALUES 

The Wigner and the Wishart ensembles belong to a class of ensembles whose joint 
probability density distribution of the eigenvalues has the form 



P(x 1 ,x 2 , 



K N exp 







N 



fc=i 
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j>i 



(13) 



where N is the number of eigenvalues, (3 is Dyson index and Kn is a normalization 
constant. In ( II 3D . V(x) is a confining potential which makes the above distribution 
normalizable by keeping the eigenvalues inside a potential. The exact expression for the 
eigenvalue density obtained integrating all variables but one is given by[3j 



N-1 



p(x) = exp[-V(x)]J2P^ 



x 



(14) 



where the P n (x) are normalized polynomials orthogonal with respect to the weight 
exp[— In fact, this is the density of the unitary ensemble, (3 = 2, but extra terms 
have to be added in the case of the orthogonal and the symplectic ensembles. The above 
sum can be reduced to the contribution of the last term using the Christoffel-Darboux 
relation 



N-1 



*£P n (x)P n (y) 



k N -i Pn{x)Pn-i{v) ~ Pn{v)Pn-i{x) 



x 



(15) 



where k n is the highest coefficient of P n (x), in the limit y — > x, (IT5|) becomes 

N-l , 

E P nW = - P^PW*)] . (16) 

The density described by Eq. ( 1141) . in the case we are interested in of Wigner's and 
Wishart's ensembles, has a central part with wiggles separated by inflection points, 
namely where the curvature changes sign, from the decaying tails at the edges. In the 
following, we use Eq. (fT6l) to show that, in the central bulk of the spectra, the density 
for these two ensembles can be written as a sum 

p{x) = p s {x) + p f (x) (17) 

of a smooth leading term, p s {x), and a fluctuating term, Pf{x), which vanishes in the 
limit of large matrices. From this decomposition of the density in smooth and fluctuat- 
ing terms, the parameters which define individual eigenvalue distributions are extracted. 
However, this procedure can not be used beyond the inflection point that is for the first 
and the last eigenvalues distributions. But, we remark that for these extreme eigenval- 
ues, the tail of their densities coincides with the density itself and does not need to be 
calculated. With this procedure, we are able to obtain the individual distribution of all 
eigenvalues. 



A. Wigner ensembles 

The joint distribution of the eigenvalues of the Gaussian random matrices for the 
orthogonal, unitary and symplectic ensembles is 

( P N \ 

P(x!, x 2 , x N ) = Z N l exp I -- E x l I II I x j - x i f, ( 18 ) 

V k=l / j>i 

where 

Z N tf) = (2 7 r)^r JV/2 ^ JV(JV - 1)/4 % + +f/2) ■ (19) 

Therefore, from Eq. ( Tl~3|) . the confining potential for this ensemble is the parabola V(x) = 
x 2 /2 and the orthogonal polynomials are the Hermite polynomials H n (x). 

1. Unitary ensemble (f3 = 2) 

The eigenvalue density is[3[ 

N-l 

P (x) = J2<f>l(x), (20) 




where 



= exp(-xV2)^(x) 
A/n!2 n v /7r 

satisfies the equation 

+ (2n + 1 - x 2 ) n (x) = 0. (22) 

which, in appropriate units, is the Schrodinger equation for the quantum harmonic oscil- 
lator. 

To be able to later calculate individual distributions of i.i.d. random variables with 
the above density, Eq. (|2"UI) . we define the counting function 



N(x) = / dx'p(x'). (23) 
Jo 

To calculate this function the recurrence relation 



4>n{x) = \ -X(j) n -l(x) ~ \ -4> n ~2{x) (24) 

v n V n 



which follows from the known recurrence relation 



H n (x) = 2xH n _x(x) - 2(n - l)H n . 2 (x) (25) 



of the Hermite polynomials can be used. Integrating the square of (1241) and using the 
relation 



(j)' n {x) = -x<f) n {x) + ^n^ix), (26) 
derived taking the derivative of f l21|) together with 

H' n = nH n _ 1} (27) 

we deduce the recursion relation 

x r 1 r x n — 1 P x 
dt<&{t) = --<t>U{x) + - / dt^it) + / dt<g_ 2 {t) (28) 

which provides an efficient way to calculate the counting function starting from the ground 
state function (f>o(x). 

Before analyzing the case of large matrices, it is instructive to consider the simple case 
of matrices of size N = 2 for which Eq. fTSOj) gives 



,,(,-) = ^^(l + 2x a ). (29) 
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The probability of the largest eigenvalue be less than a value x is obtained by integrating 
the joint distribution P{x\,x_) in the two variables X\ and x 2 from — oo to x. Taking 
then the derivative of this probability, we find that the density distribution of the largest 
eigenvalue is 



;i + 2x 2 )(l-erf(x))-— GXP( ^ 



7T 



(30) 



For the smallest eigenvalue distribution it is simpler to determine it by taking the differ- 
ence p{x) — F(0,x). For the uncorrelated case (see previous section), with n — 0, 1 and 
p(x) given by ([29|) the density distributions of the largest and the smallest are given by 



Fu(n,x) 



p(x) 



1 + (-l) n erf(x) - 



xexp( 



-x 



71 



(31) 



In Fig. 3, these distributions are shown together with the density. One can clearly see 
the reduction produced by the correlations in the range of fluctuations of the eigenvalues. 

Turning now to the case of large matrices, by using the ( TTBj) . (T2"7]) and (f2"T|) we rewrite 
flUD as 



(32) 



p( x ) = V T [0jv(»<^-iO) - 0^0)0^-1 (>)] . 



If the harmonic oscillator functions are expressed in terms of amplitude and phase as 



4> n (x) = A n (x) cos8 n (x) 



(33) 



then, when substituted in ( |32|) . the density becomes a sum of cosines and sines whose 
arguments are either addition or subtraction of the function phases. We observe that as 
the function <p n (x) has n zeros its phase 6 n (x) must have a range of variation of order 
mi. Therefore, by subtracting two adjacent phases we get a function whose variation is 
smaller than tt producing therefore no wiggles. We can then argue that the smooth part 
of the density, p s (%), is obtained by collecting the terms in which the phases subtract 
while, the fluctuating part, Pf{x), comes from those in which they add. Explicitly, this 
leads to the expressions 
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(35) 

(all quantities 6 n 's and A n 's above and others below are x-dependent, but to keep the 
notation less heavy, this dependence is often dropped). An asymptotic expression for 
these two density terms can be deduced from the semi-classical formalism described in 
appendix A. We find that, asymptotically, the harmonic oscillator function is given by 



<Pn{x) 



where 



2 COS [in{x) — n/2] 7T 

?r (2n + 1 - x 2 ) 1/4 
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2n + 1 



(37) 



is the classical mechanical action obtained integrating the momentum p = \/2n + 1 — x 2 . 
In ( 136]) . the phase nn/ 2 has been introduced to fix the parity of the function. Eq. f )36|) 
together with ( )37|) determine the phases and amplitudes in ( )34|) and ( |35l) Substituting 
them and neglecting the derivatives of the amplitude, we find, after neglecting higher 
order terms, that the density takes the simple analytic form 



p(x) = p w {x) 



/2iVcos[iV-2£(x)]7r 



where 



2tt 3 p 2 
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(38) 



Pw{x) 

is Wigner's semi-circle law and 



x z 




(39) 
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arcsm 




(40) 



In deriving these equations, N was assumed to be large enough to treat indexed quantities 
as continuous functions of the indices in such a way that the approximations f(N) + f(N— 
1) = 2f(N - 1/2) and f(N) - f(N - 1) = f'(N - 1/2) can be made. 

Eq. ( l38l) shows that averaging out the wiggles the fluctuating term vanishes and the 
density becomes the semi-circle. The quantity £(x) is the average stair-case function from 
which the so-called unfolded spectrum is calculated. This transformation leads to a new 
spectrum with density 



p(0 



Pw(x) 



/2iVcos(A^ - 2f)7r 



(41) 



whose average is equal to one. Since £(x) is a counting function, in this variable the 
wiggles are equally spaced with an unit distance between them. In Fig. 4, the density 
calculated using the approximated expression Eq. ( 138]) is compared with the density 
calculated with Eq. (!20l) . One can see that with the exception of the wiggles at the very 
edges, those at the bulk are well described by the asymptotic density. 

Eq. ( 14 ip gives analytical expressions for the smooth and the fluctuating parts of the 
density which Fig. 4 shows are very precise at the bulk of the spectrum. To extend this 
precision up to the edge, removing the singularities at the classical turning points, an 
improvement is needed to go beyond the asymptotic expression Eq. (|36|) . To this end, 
one must also consider the exact second independent solution of the equation expressed 
as 



[x 



AJx) sinO J x) 



(42) 



From the pair of independent solutions, modulus and phase of the functions can be 
extracted. This second independent solution can be determined by integrating from the 
origin, x = 0, the differential equation with initial conditions provided by the Wronskian 
relation 



W((j) n , 4> n ) = <j> n {x)4/ n {x) - (f/ n (x)(f> n (x) = 2/tt. 



(43) 



In fact, since solutions of the equation can be constructed with a definite parity, for even n, 
we take (f> n (x) to be odd such that n (O) = and, from (jSP , <fi' n (0) = 2/ii(f) n (0). Inversely, 
for n odd, 4> n {x) is taken to be even, 0'„(O) = and, from (03]), n (O) = — 2/7r0' n (O) . 
Once the pair of solutions is determined, modulus and phase are obtained as 
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and their derivatives as 
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The exact expression of the scaled density becomes 



p(x) _ ^ B cos(9 N + 6>at_i + 



where 
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and 



a , A N A N _x(9' N - v N _u 

9 = arctan — — . (50) 

Guided by intuition and numerics, let us make the ansatz that the eigenvalue density 
in the scaled variable can be decomposed as 



9' N _ } , 



k=l ^ 27la 



where, according to Eq. (138 
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with k = 1, 2, N. Since the quantities v and a may depend on position, each term in the 
above sum is not a true Gaussian, however they can be considered as nearly Gaussian 
distributions as that dependence is expected to be weak. In order to determine this 
dependence, we turn the summation in (15"TI) into an infinite sum and rewrite it as 
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where 
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This R quantity is expected to affect only the tail of the distribution of the extreme 
eigenvalues which, as explained in the beginning of this section, practically coincide with 
the total density. Therefore, this reminder can be neglected even for relatively small 
values of N. The first term in the right hand side of (153]) can be transformed into an 
infinite sum of integrals using the Poisson sum formula 



S2 f(t + n)= ^2 F(27rm)exp(27rimt), 

n=— oo m=— oo 

where F(s) is the Fourier transform 
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Then, after performing all integrals we obtain 
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p{v) = l + 2^(-l) m exp [-2(imia) 2 ] cos 
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(57) 



where the term R in (153]) was neglected. Due the presence of the exponential factor, the 
sum is dominated by its first term m = 1. Comparing this term with the oscillating term 
in Eq. (I4ip . we find that the phase v and the variance a depend on their positions as 
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At the bulk of the spectrum, the phases take their asymptotic values 9 n = (£„ — nj < £) r K 
and 9 = n/2 such that u(x) = £(x) and 
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Once the variances have been determined, the decomposed density in the actual spectrum 
variable is 



N 1 
p(x) =p s (x)J2 7k 2 , - s 



exp 



{v(x) - v k y 
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In Fig. 5, these N individual distributions are compared with the distributions ob- 
tained by performing numerical simulations for matrices of size N = 20 with a very good 
agreement. In Fig. 6, these density distributions are compared with those of uncorrelated 
variables with the same density exhibiting the great effect of the correlations. 

Motivated by the good agreement between simulations and nearly Gaussian distribu- 
tions, we compare the nearly Gaussian distribution at the edge with Tracy and Widom's 
prediction for the largest eigenvalue. They proved that, when N — > oo, in a new variable 
s defined by the linear relation 

x = v/^V+^t^-^, (62) 



2 i/2jyi/6' 

the distribution probability, E 2 (s), of the largest eigenvalue of the unitary ensemble, 
— 2, is given by [gj 



exp 



(x — s)q 2 (x)dx 



where q(s) satisfies the Painleve II equation 

sq + 2q 3 



(63) 



(64) 



with boundary condition 



q(s) ~ Ai(s) when s — > oo, 



(65) 



where Ai(s) is the Airy function. 

In Fig. 7, the distribution of the largest eigenvalue of matrices of size N = 20 obtained 
performing numerical simulation is compared with both: Tracy- Widom density distribu- 
tion, E' 2 (s) and our nearly Gaussian distribution. Both give a reasonable fit although 
N = 20 can be considered a relatively small size. In Table 1, the cumulants of the two 
distributions are shown (l7| and, as one would expect, the figures point that ours is indeed 
more normal. 





mean 


variance 


skewness 


kurtosis 


Tracy- Widom 


-1.77109 


0.9018 


0.224 


0.093 


nearly Gaussian 


-1.829 


0.9066 


0.114 


0.074 



Table 1: Cumulants for the unitary case 



2. Orthogonal ensemble (13 = 1) 
For the orthogonal case, we have to add to (|T4|) the term 



70) 



N f°° 1 

— (j) N -i{x) I dt-sgn (x - t)(j) N (t) 



<Pn-i{x) / dt(f) N (t), (66) 



where sgn (x) is the sign function, if iV is even, and a further term 



oc 



Mx)/ / <f> N {t)dt, (67) 



if A" is odd@. 

As done in the unitary case, it is instructive to start discussing the decomposition of 
the spectral density of matrices of size N = 2. This will illustrate the differences between 
the unitary and the orthogonal cases. Evaluating the integral in (I66I) we find that the 
extra term is 

, . 2x 2 exp(-x 2 ) exp(-x 2 /2)x f x , . t 2 . 
7 (x) = ^ -+ 1 ' ' / dtexp{--). (68) 

Adding this term to (I2"9~j) . its first term cancels the second term in (I29p and the density 
becomes 

. , exp(— x 2 ) exp(— x 2 12) „. x . , . 

P W = ^^ + ^ r ^erf(^). (69) 

The important point is that this cancellation of terms removes the wiggles in the uni- 
tary density in such a way that the orthogonal becomes a flat function. The individual 
distributions of the two eigenvalues are easily calculated to be given by 

m x) = V?5E<ZE2 erfc(-^) + 5^ (70) 

for the greater while distribution of the smaller is obtained by subtracting the above 
F(0,x) from Eq. flo^|) . For the uncorrelated case, with n = 0, 1 and p(x) given by 
the density distributions are 



F u (n 1 x) = ^Y 



1 + (-1)" ( erf(x) - eXp( ^ 2/2) erf(^; 



(71) 



Fig. 8 shows the spectral density of matrices of size A^ = 2 of the orthogonal ensemble 
decomposed in eigenvalue individual contributions. In contrast to the unitary case , Fig. 
3, in the orthogonal case, the two individual contributions add in such a way that density 
becomes flat at the top. 

Turning now to the case of matrices of large sizes, we have to calculate the integral 

I N (x)= [ dt(j) N (t). (72) 
Jo 

Starting with the derivation of an asymptotic expression for it, we use the differential 
equation satisfied by the 4>n{x) functions to rewrite the it as 

We recall that the denominator in the above integrand is the square of the classical 
momentum supposed to be large. Therefore integration by parts can be used to obtain a 
series in inverse powers of the momentum or, equivalently, of the density, whose the first 
two terms are 



/iY(x) -"2iV+l-^ + (2iV + l-xT- ( } 

Substituting this term in Eq. ( )66|) and replacing the functions by their semi-classical 
approximation, we find that the first additional extra term cancels the oscillating term 
of the unitary case. This canceling makes necessary to take into account higher order 
terms, which can be done using the expansion 

(2JV ± 1 - x 2 Y = (2N - x 2 Y ± fi(2N - x 2 y~ l . (75) 
By doing this, we end up with the following asymptotic expression for the density 



p{x) = pw -2^ + 2^; cos (5 " yj « + -s^ sm - -2) 71 (76) 

in the orthogonal case. The first two terms in (I76p correspond to the smooth part of the 
density while the two last ones to its fluctuating part. 

In order to go beyond this asymptotic expression, removing, as done in the unitary 
case, its singularities, we assume, that the integral In(x) can be written as 

In(x) = Q\(x) cos9 N (x) + Q 2 (x) sin9 N (x), (77) 

where 9n{x) is the phase of the function 4>n(x) and Qi(x) and Q 2 (x) are smooth functions 
of the position. These functions can be determined considering the function related 
to the second independent solution as I' N (x) = 4>n(x). From the asymptotic analysis, we 
deduce that it satisfies, at the origin, the condition 



I N (Q) = - . f 7 i 

K ' 2N + 1 v 



and can be written as 



In(%) — Qi(x) sm9 N (x) — Q2(x) cos9 N (x). (79) 
Eqs. (177)) and (179"]) can be inverted to give 

Qi(x) = In(x) cos9n(x) + In{x) sin^Ar(x), (80) 

and 

Q2{x) — In(x) sin 9n(x) — In(x) cos9n(x). (81) 

Once these two functions are determined the integral Eq. (I77|) is expressed in terms of 
amplitudes and phase and can be replaced in the additional term j(x) which can also be 
decomposed as a sum of smooth and fluctuating terms as 

[W 

7s = \ —A N ^ [Q 1 cos{9 N - e N - X ) + Q 2 sin{9 N - 9 N ^)] (82) 
V o 

and 




7/ = \l — a n-i [Qi cos(8 N + 9 N -x) + Q 2 sm(9 N + 9 N ^)] . 
By adding 7 S and 7/ to Eqs. flMj) and fl35|) respectively the scaled density is 

Bcos(9 N + 9 N -! + 9) 



P 



pl s (x) 



Pis 



where 



Pi s (x) = p s (a;) +7 s (x), 



(83) 



M) 



(85) 



,'JV 

B = \ — A N A 




g 2 



and 



= arctan 



9' N — 9'r 



N-l 



Q2/A 



N 



A' N /A N - A'^JAn-i + Qi/A N 

The formalism developed in the unitary case still applies and, with the above 
sions, the scaled variable and the variance are deduced to be given by 



(86) 

(87) 
expres- 



u(x) 



7/V 



+ 0JV-1 + 



7T 



and 



a 2 (x) 



2tt 2 



2tt 



In 



N 
2" 



2pi s (x) 



At the bulk, from Eq. ( ITS"]) , the scaled variable becomes 

t'(x) = + arctan — 

_4:irpw{X 

and the density smooth term and the amplitude are given by 



Pis Or) = Pw{x) 



1 



2n 2 p w {x) 



and 



(89) 



(90) 



(91) 



B(x) 



V2N 



1 + 



9x 2 



I6ir 2 p 2 w 



(92) 



respectively, which replaced in Eq. f|89|) gives the asymptotic expression of the variance. 

In Fig. 9, the density distributions of individual eigenvalues of matrices of size N = 20 
are shown together with distributions obtained performing numerical simulations. The 
result is good specially at the bulk of the spectrum. 

For the largest eigenvalue, Tracy and Widom predict that for the orthogonal ensemble 
((3 = 1), the probability distribution function in the same scaled variable s of the unitary 
case is given by 



[£!(*)]* ^OOexp i-n(s)] 



(93) 



where 

;>oo 

fjL(s) = / q(x)dx. (94) 

J s 

In Fig. 10, this prediction is compared with our nearly Gaussian density distribution 
and with results of simulations. It is clear that the nearly Gaussian distribution give a 
better description. However, comparing the cumulants|l7j] depicted in Table 2, we see 
that actually, apart from a shift to the right of the Tracy- Widom distribution, the two 
distributions are quite alike. 





mean 


variance 


skewness 


kurtosis 


Tracy- Widom 


-1.20653 


1.2580 


0.293 


0.165 


nearly Gaussian 


-1.382 


1.264 


0.325 


0.067 



Table 2: Cumulants for the orthogonal case 



B. Wishart matrices 

Consider a rectangular matrix X of size (MxN) whose elements are sorted indepen- 
dently from a Gaussian distribution, a Wishart square matrix W of size N is then defined 
by taking the product W = X^X. It can be shown that for M > N, the joint probability 
distribution of the positive eigenvalues of the random matrices W, for the three symmetry 
classes, are given by 

P(X 1: X 2 , ...,X N ) = jfjyexp(-^y^37 fc ) Y\x? + Y[\ X 3~ X if- ( 95 ) 

k=l i=l j>i 

From (TT3~]) . the positive eigenvalues of this ensemble are confined by the potential V(x) = 
x—(l+M—N—2/ 0) log(x) and the polynomials are the generalized Laguerre polynomials. 
For the unitary case, (3 = 2, the eigenvalue density is 

N-l 


where, with a = M — N, 



r n {x) = J-^—-exp(-x/2)x a / 2 L a n (x) (97) 



in which the L%(x) are the Laguerre polynomials 



^w-B-y (n ^ )!jr «») 



From (jgg|) . we derive 



dx 



T a+1 



which used in (fl6j) yields 

p(x) = exp(— x)x 



T{N) 



[L%_ 1 (x)L'^ n 1 (x) — L%(x)L'^ h _} 2 (x)j . 



(99) 



(100) 



T(N + a) 

Inverting (l97l) to express the polynomials in terms of the ^-functions, the density becomes 



a+l 
N-2 



X 



(101) 



Finally the derivative of Eq. ( 1971) gives the relation 



dx 



a — x 
2x 



T Tl 



which used in (110 II) allows to write the density as 



d 



d 



(102) 



(103) 



an expression now ready for our analysis. We start by remarking that the ip n (x), in 
appropriate units, are the functions which satisfy the Schrodinger wave equation of the 
hydrogen atom 



d 2 i\) dip 
dx 2 dx 



1 

U+ 2 



(a — x) 
Ax 



ip = 



(104) 



(from now on, the superscript a will be omitted). As in the Wigner case, we write these 
functions in terms of amplitude and phase as 



ip n (x) = A n (x) cos 9 n (x) (105) 

and substitute, in the expression of the density, the two functions to extract its smooth 
part 



' N(N + a) A 

n, = \i — ^ — -AnAn^ 



and its fluctuating part 



A 



N-l 



A 



N-l 



A' 

A N 



COS(U N 



- e 



N-l. 



+ 



N 



9' 



N-l 



sin! V N 



- 9 



N-l i 



(106) 



Pf 



N(N + a] 



A N A 



N-n-N-l 



A 



N-l 



A 



N-l 



A' 

A N 

A N 



cos(9 N + 9 N - X ) + {9' N - 9'^) sin{9 N + 9 N _ t ] 

(107) 

In order to apply the semi-classical formalism, we first transform equation (II 04ft in 
the differential equation 



f{ -f*> 4- - x2 + 2(2 " + °+ 1)l - a2 + 1 (^) = 0. (108) 
satisfied by the function y/xip(x). In this equation, the associated classical moment is 

p(x) = ^\f(x-x 1 )(x 2 -x), (109) 

where 

Xl>2 = 2n + a + 1 =f y/(2n + a + l) 2 + 1 - a 2 . (110) 
and the asymptotic wave function is 

/ ( \ /2 cos(£ ra - tt/4) 
^ = Vh^-x^x-x^ 



where 



1 



-V^arctan yfg^ + (x 2 + a*) a rccos (^g^) 

+2a/ (x 2 - x)(x - xi) 



(112) 



Substituting this approximate wave function in the smooth and the fluctuating parts of 
the density and neglecting derivatives of the amplitudes, we arrive at the asymptotic 
expression 

p{x) = Pmp(x) - 167 ^2 p 2^ (x ) cos[2£(x)]tt, (113) 
where the first term is the Marchenko-Pastur density [3] 

Pmp(x) = ^—\ / (x + - x)(x - x_), (114) 
lixx 



K 

N 

cosine argument is the counting number function 



in which, with c = yjj, x± = iV(c ± l) 2 . Similarly, the function appearing in the 



— A./x_x + arctan \j x , x \ + (x+ + a;_) arccos ( — — ) 

V + V x-{x+-x) 1 V + 1 ) V x+-x- J (115) 

+2a/ (x + — x)(x — x_) 

associated to the Marchenko-Pastur density. As in previous case, approximations were 
made by treating indexed quantities as continuous functions of their indices. Eq. dl 13[) 
shows that the Marchenko-Pastur density plays, for the Wishart matrices, the same 
role the Wigner's semi-circle does for the Gaussian ensembles and, by averaging out the 
wiggles produced by the oscillating term, it is obtained. We remark that the above 
expression for the next to leading order term of the asymptotic, has an analog structure 
to that of the Gaussian expression Eq. ( I4"T1) . namely, both are, basically, the ratio between 
the superior limit and the cubic of the asymptotic density. As before, we use the counting 
function £ as independent variable with density 



In Fig. 11, we compare the above density approximation with its exact expression. It is 
seen that it is indeed a very good approximation with the exception of the regions near 
the two borders. 

As in the Wigner case, a complete decomposition of the density which includes eigen- 
values at the border can be achieved using the second independent solution of the wave 
equation which can be obtained from the integral representation 



2«(_i)t 

a/2 



1XX 



n + a)! 



ft! 



d9 cos" 1 9 cos 



xtan^ 



(2ft + a + 1)9 



(117) 



of the first solution. That this function, with a > 1, is equivalent to (1971) can be seen 
by changing the integration variable to u — tan 9 that transforms the integral in ([TTTJ in 
the complex integral 



2 a (-l) n (n + a)\ 



2irx a / 2 



rv. 



-0 



2n+a+l 



(118) 



(ft - 

which performed by residues reproduces (197fl . This integral representation suggests that 
the other independent solution is provided by replacing in the integrand the cosine of the 
oscillating factor by the sine. However, this is not enough, and yet another term has to 
be subtracted to construct the solution 



2«(-l)» l(n + a)\ f lo d9cos a - 1 9sm[^-(2n + a + l)9] + 



TlX 



a/2 



111 



J™dB cosh 1 *- 1 9 exp 



;ii9) 



tanh 6 



(2ft + a + 1 



Now, we have a pair of independent solutions and modulus and phase are determined 
through the relations 



and 



with derivatives given by 



A, 



9 n = arctan(-0„/'0n) 



(120) 



(121) 



A n A' n = ip n i/j' n + Tp n ip' n 



and, using the Wronskian W(ipip)(x) = 



n _ n ' 
7T 



;i22) 



(123) 



In calculating these quantities the two recurrence relations 



Zn-2 — 



y/(n- l)(n + a- 1) 



(2n + a - 1 - x)Z n _i - yjn{n + a)Z n 



(124) 



and 



;r- 



dZ n f 2n + a — x 
dx 



/ 2n + a — x\ ^ n i—. - 

f g ) ^ " V^ + ^J^n-l, 



(125) 



where Z n denotes any one of the two solutions are useful. 

The precise expression for the decomposition of the density in smooth and fluctuating 
parts is 



p{x) 

Ps{x) 



1 - 



Bcos(9 N + 9 N _ l -9) 



(126) 



where 



B 



and 



9 = arctan 



A N A 



N-rt-N-l\?N ~ Q'n-1 



(127) 



(128) 



. An-iA' n — AnA' n _ 1 

Turning now to the decomposition of Wishart eigenvalue density in a sum of individual 
nearly Gaussian distributions, we write 



JV-l 



: exp 



2a 2 



where 



v k = - + k. 



(129) 



(130) 



with k = 0, 1, 2, N — 1. As was done in the Wigner case, we turn the above sum into 
an infinite sum which can be transformed, after neglecting the rest, into the infinite sum 



p(u) = 1 + 2 ^(-l) m exp [-2{-Kmaf] cos [2nm(u)] , 



131) 



m=l 



using Poisson sum formula. The dependence of the phase and the variance on their 
positions along the spectrum is then determined to be given by 



v{x) = 



On + 9_n-i — 9 
2^ 



(132) 



and 



a\x) 



2tt 2 



In 



' B(x) ' 



(133) 



At the bulk, these quantities approach their asymptotic expressions u(x) = £(x) and 



a 2 (x) 



2tt 2 



In 



2n(2x) 2 ^p MP (x) 



x 



1/3 



(134) 



In Fig. 12, the individual eigenvalue distributions of N = 20 Wishart matrices are 
shown. It is seen that also for this ensemble our formalism based on nearly Gaussian 
distributions give account of the results obtained performing numerical simulations. 

For the largest eigenvalue the prediction is that its distribution follows, for large N, 
the Tracy- Widom distribution in the scaled variable 0] 



x — X- 



x 



2./3. 



(MN)- 



-1/6 



(135) 



In Fig. 13, the two theoretical predictions are compared with results from numerical 
simulations. Finally, in Fig. 14, the density distributions of the smallest eigenvalue is 
magnified to show that the nearly Gaussian distribution gives a good fit to the numerical 
simulation. 



IV. CONCLUDING REMARKS 

We have performed a decomposition in individual contributions of the exact density 
of eigenvalues of matrices of the unitary and the orthogonal Gaussian ensembles and of 
the unitary matrices of the Wishart ensemble. This decomposition works for relatively 
small values of matrix sizes and shows that eigenvalue are well described by nearly Gaus- 
sian distributions. As the matrices sizes increase, these distributions tend, at the bulk of 
the spectra, to true Gaussian such that the exact results of Refs. jsl, |To[ are recovered. 
These results should apply to spectra of systems whose classical analog are chaotic and 
be experimentally tested. The present analysis should also work for ensembles connected 
to others orthogonal polynomials. Recently, invariant non-ergodic ensembles have been 
introduced whose ensemble densities are averages over the Wigner's semi-circle 19[ , in 
the Gaussian case, and the Marchenko-Pastur density (2oT]. in the Wishart case. In the 
Gaussian case, the effect of non-ergodicity on the Tracy- Widom distribution has been 



investigated [21| . It seems possible to extend this study to the individual distributions of 
the eigenvalues at the spectral bulk. This is interesting in connection with the behavior 
of the eigenvalues of the embedded ensembles^, 22[. On the other hand, the high de- 



velopment the experimental physics achieved in the field of cold atoms, should perhaps 
reach the point in which individual atoms can be observed. In this case, the behavior 
of each atom of a Girardeau gas can be measured and the decomposition in individual 
contributions be checked. 

This work is supported by a project CAPES-COFECUB and by the Brazilian agencies 
CNPq and FAPESP. 



Appendix A: Semi-classical wave function approximation 



Consider a wave function \l/(x) that satisfies the wave equation 



d 2 m 

dx 2 



+ p 2 (x)^ = 0, 



(Al) 



where the momentum p(x) is such that p 2 (x) > in the interval £2). Let us search a 
solution of PAID in this interval of the form 



V(x) = M(x) cos 0(x). (A2) 

Substituting ( 1A2|) in flAlj) and neglecting the second derivative of the modulus, we find 
that modulus and phase satisfy the equations 



(A3) 



and 



M(3" + 2M'(3' = (A4) 

with solutions f3 = f* dxp(x) + (5q and M = Cp~ x l 2 . Finally the constant C is fixed by 
the normalization condition J ^> 2 (x)dx = 1 that gives 



C 



1 f X2 dx[l + cos 2/3 (x)] 



x-l 



p(x) 



- -1/2 




r 2 dx ' 


-1/2 




2 







(A5) 



where, consistently with the semi-classical method, the oscillating term in the integrand 
is neglected. 

As a second order differential equation, Eq. AA1H has a second independent solution 
which in the same approximation we assumed to be given by 



ty(x) = M(x) sin/3(x). 
These two independent solutions satisfy the Wronskian relation 

in the same order of approximation. 
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Figure Captions 

Fig. 1 Distributions of A" = 20 individual Gaussian variables with density 
N exp(-x 2 ) / y/n . 

Fig. 2 For A" = 100, histogram of the largest Gaussian variable scaled as t(x) = 
yerf(x) is plotted together with the predicted Gumbel density distribution. 

Fig. 3 Decomposition of the N = 2 density of the unitary ensemble (the black line) 
in individual contributions: eigenvalues (blue lines) and uncorrelated random variables 
(red line). 

Fig. 4 The relative difference (p aS ym ~ p)/p between the exact density p and its 
asymptotic is plotted versus the unfolded variable £. 

Fig. 5 For the unitary case, the A" individual distributions calculated with our formal- 
ism complete (blue line) and asymptotic (red line) compared with the result of numerical 
simulations (black line) for A" = 20. 

Fig. 6 Decomposition of the A^ = 20 density of the unitary ensemble (the black line) 
in individual contributions: eigenvalues (blue lines) and uncorrelated random variables 
(red line). 

Fig. 7 Theoretical distributions of the largest of A^ = 20 eigenvalues compared with 
numerical simulation (blue line): nearly Gaussian (black line) and Tracy- Widom (blue 
line) and asymptotic (red line) compared with the result of numerical simulations (black 
line) . 

Fig. 8 Decomposition of the N = 2 density of the orthogonal ensemble (the black line) 
in individual contributions: eigenvalues (blue lines) and uncorrelated random variables 
(red line). 

Fig. 9 For the orthogonal ensembles, A^ = 20 individual distributions calculated with 
our formalism complete (blue line) and asymptotic (red line) compared with the result 
of numerical simulations (black line). 

Fig. 10 For the orthogonal ensemble, largest eigenvalue density distribution for A" = 
20: nearly Gaussian (black line), Tracy- Widom (red line) and simulation (blue line). 

Fig. 11 The relative difference (p aS ym — p)/p between the exact density p and its 
asymptotic is plotted versus the unfolded variable £. 

Fig. 12 Decomposition of the A^ = 20 density of the unitary ensemble (the black line) 



in individual contributions: eigenvalues (blue lines) and uncorrelated random variables 
(red line). 

Fig. 13 For the Wishart ensemble, largest eigenvalue density distribution for N = 20: 
nearly Gaussian (black line), Tracy- Widom (red line) and simulation (blue line). 

Fig. 14 For the Wishart ensemble, the smallest eigenvalue density distribution: nearly 
Gaussian (black line) and simulation (red line). 



